First-order phase transition in the tethered surface model on a sphere 
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We show that the tethered surface model of Helfrich and Polyakov-Kleinert undergoes a first-order 
phase transition separating the smooth phase from the crumpled one. The model is investigated 
by the canonical Monte Carlo simulations on spherical and fixed connectivity surfaces of size up 
to N = 15212. The first-order transition is observed when N > 7000, which is larger than those in 
previous numerical studies, and a continuous transition can also be observed on small-sized surfaces. 
Our results are, therefore, consistent with those obtained in previous studies on the phase structure 
of the model. 
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I. INTRODUCTION 

Considerable progress has been made in understand- 
ing the phase structure of the model of Helfrich and 
Polyakov-Kleinert d 13, 0, S IE IE 0, S, S EO, E3, E3 ; we 
will abbreviate this to the HPK model. The HPK model 
describes not only two-dimensional surfaces swept out by 
a rigid string |7| but also biological membranes such as 
human red-blood cells and artificial vesicles d 0, H El ■ 

It is now widely accepted that the tethered model 
of HPK undergoes a continuous phase transition. The 
Hamiltonian of the HPK model is given by a linear com- 
bination of the Gaussian term S\ and the bending energy 
term S2 : S = S\ + bS 2 , where b is the bending rigidity. 
52 is ordinarily defined by using the unit normal vec- 
tors for the triangles. The large-Z? expansion E3 pre- 
dicts that the HPK model undergoes a finite-6 contin- 
uous transition between the smooth phase in the limit 
b —> 00 and the crumpled phase in the limit b — > 0. Nu- 
merical studies performed so far have also focused on the 
phase transition in the tethered surface models of HPK 
[11 El E H El El and indicate that the model ex- 
hibits the continuous transition. 

On the contrary, we can also think that the model 
has a discontinuous transition. It was predicted by mean 
field analysis that the model undergoes a first-order phase 
transition [l^. I n recent numerical simulations on fixed 
connectivity surfaces [H |2l| , it was also suggested that 
the phase transition is of the first order. The Hamiltonian 
of the model in |2fJ includes the Lennard-Jones potential 
serving as the Gaussian energy for the HPK model. The 
bending energy in |2l| is very similar to the one for the 
HPK model. These numerical results, therefore, strongly 
suggest that the HPK model undergoes a discontinuous 
transition. 

However, little attention has been given to whether 
a discontinuous transition is observed in the tethered 
surface model of HPK, whose Hamiltonian includes the 
bending energy of the form 1 — • xi 3 ■ , where is the 



unit normal vector for the triangle i. We will call this 
form of energy as the ordinary bending energy from now 
on. Therefore, in order to confirm that the phase transi- 
tion of the HPK model is of the first order, we need to 
study further the tethered surface model defined by the 
ordinary bending energy. 

In this paper, we numerically study the tethered model 
on a sphere, whose Hamiltonian is given by S=Si+bS2, 
where Si is the Gaussian energy and S2 is the ordi- 
nary bending energy described above. This Hamilto- 
nian has been widely accepted and investigated as a 
discrete model of HPK. Although a fluid surface model 
[H HI HI HI E3 defined on dynamically triangu- 
lated surfaces is very interesting and should be investi- 
gated further on larger surfaces, we will concentrate on 
the fixed connectivity tethered model in this paper. 

We will show the first numerical evidence that the or- 
dinary tethered model undergoes a discontinuous transi- 
tion on a sphere. The gap of the bending energy is clearly 
seen on surfaces of N > 7000, and cannot be seen on the 
smaller surfaces. It must be emphasized that the results 
are not contradictory to previous ones, as the continu- 
ous transition can also be observed in our simulations on 
smaller surfaces. 



II. THE MODEL AND MONTE CARLO 
TECHNIQUE 

The partition function of the model is defined by 



N 

Z = J J[dXiexp[-S(X)], 

i—1 



S(X) =St + bS 2 



(1) 
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where b is the bending rigidity, N the total number of 
vertices. The center of the surface is fixed to remove the 
translational zero mode. The self-avoiding property of 
surfaces is not assumed in the integrations dXi in R 3 . 
The symbols Si , S 2 in Eq. (JTJ denote the Gaussian en- 
ergy and the bending energy, which were already intro- 
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duced in the Introduction and are defined by 



Si 



Xj) , S2 



*i) , (2) 



where JfJ^ denotes the sum over bonds (i, j), and ^ ^ 
the sum over triangles i,j sharing a common bond. The 
symbol in Eq. J2J denotes the unit normal vector of 
the triangle i, as was introduced in the Introduction. 

The canonical Monte Carlo (MC) technique is used 
to update the variables X . The new position X[ of the 
vertex i is given by X[ — Xi + AX , where AX is chosen 
randomly in a small sphere. The radius of the small 
sphere is chosen at the start of the MC simulations to 
maintain about a 50% acceptance rate. The new position 
X[ is accepted with the probability Min[l, exp (— AS)], 
where AS is given by Z\S' = S'(new) — 5(old). 

The minimum bond length is not assumed. On the 
contrary, the minimum area of triangle is assumed to be 
1CP 6 x Aq, where Aq is the mean area of triangles com- 
puted at every 250 Monte Carlo sweeps (MCS) and is 
almost constant throughout the MC simulations. The 
area of almost all triangles generated in the MC simula- 
tions is larger than the lower bound 10 -6 x Aq. 

We use a random number called Mersenne Twister 
in the MC simulations. Two sequences of random num- 
ber are used; one for a 3-dimensional move of vertices 
X and the other for the Metropolis accept /reject for the 
update of X. 

The surfaces, on which the Hamiltonian in Eq. is 
defined, are obtained by dividing the icosahedron, and 
hence, are uniform in the co-ordination number. By 
dividing every edge of the icosahedron into L pieces of 
the same length, we have a triangulated surface of size 
N = 10L 2 + 2. These surfaces are thus characterized by 
-/V5 = 12 and N@ = N— 12, where N q is the total number 
of vertices with co-ordination number q. Hence we have 
surfaces in which 12 vertices are of qi = 5, and all other 
vertices = 6. Hence, the surfaces are made uniquely 
in contrast to the Voronoi lattices constructed by using 
random numbers. 

We comment on the unit of physical quantities. The 
scale of length in the model can be chosen arbitrarily 
because of the scale invariant property of the partition 
function in Eq. l|T]l. Then, by letting a be a length unit 
(the mean bond length for example) in the model, we can 
express all quantities with unit of length by a, which is 
assumed to be a = 1 in the model. Hence, the unit of Si 
is a 2 . Let A be the surface tension coefficient, which is 
assumed to be A = 1, S in Eq. can also be written 
as S = XSi + bS2- Thus, the unit of A can be written as 
kT/a 2 , where k is the Boltzmann constant and T is the 
temperature. The unit of 6 is then expressed by kT . 
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FIG. 1: (a) The mean square size X 2 against the bending 
rigidity 6, and (b) the bending energy S2/NB against b, where 
Nb is the total number of bonds. Both X 2 and S2/NB discon- 
tinuously change at some intermediate 6. This discontinuity 
represents some discontinuous phase transition. The unit of 
X 2 is a 2 , where a is the length unit in the model, which is 
chosen to be a = l. 



III. RESULTS 

We firstly show in Fig. da) the mean square size X 2 
defined by 



where denotes the sum over vertices i. Discontinuous 
changes of X 2 are visible at intermediate bending rigidity 
b in Fig. ^a), and suggest that there is a discontinuous 
transition between the smooth and the crumpled phases. 

Figure ^b) shows the bending energy S2/NB against 
b, where Ng is the total number of bonds. We find a 
discontinuity in Si/Nb at b where X 2 discontinuously 
changes. This can be viewed as a signal of a discontinuous 
transition. 

Total number of MCS at the transition point was 5 x 
10 8 , 8xl0 8 , and 8xl0 8 for surfaces of N = 7292, N= 10242, 
and N= 15212, respectively. A relatively small number 
of MCS were done at b far distant from the transition 
point on each surface. 

Figures|2[a), (b), and (c) represent the variation of X 2 
against MCS, which were obtained at 6 = 0.77, 6 = 0.772, 
and 6 = 0.773 on the surface of size 7V = 7292. We can find 
from Fig. HJb) that there are at least two states which 
differ in size; one of them is characterized by X 2 ~ 30 and 
the other by X 2 ~ 70, which correspond to a crumpled 
and a smooth state respectively at the transition point 
of the TV = 7292 surface. The size of the surface appears 
to be stable, and this stability of size reflects a phase 
transition. Thus, we understand that the surface remains 
in the crumpled (smooth) phase at 6 = 0.77 (6 = 0.773), 
and that the transition point is close to 6 = 0.772 on the 
TV = 7292 surface. Figures Eld), (e) and (f) are those 
obtained at points 6 = 0.769, 6 = 0.77, and 6 = 0.771 
respectively on the surface of N = 10242. We can see 
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FIG. 2: The variation of X 2 against MCS obtained on the N = 
7292 surface at (a) 6 = 0.77, (b) 6 = 0.772, and (c) 6 = 0.773, 
those obtained on the N = 10242 surface at (d) 6 = 0.769, 
(e) 6 = 0.77, and (f) 6 = 0.771, and those obtained on the 
N = 15212 surface at (g) 6 = 0.768, (h) 6 = 0.769, and (i) 
6 = 0.77. 
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FIG. 3: (a) The normalized histogram h(S2) of Sz/Nb, and 
(b) the variation of S 2 /N B against MCS, on the N = 7292 
surface at 6 = 0.772, (c) h(S2) and (d) the variation of S2/NB 
on the N = 10242 surface at 6 = 0.77, (e) h(S 2 ) and (f) the 
variation of S 2 /N B on the N = 15212 surface at 6 = 0.769. 
Two distinct peaks on each h(S 2 ) indicate that the model 
undergoes a discontinuous transition. 



that the transition point for the N — 10242 surface is 
close to b — 0.77, where two states which differ in size 
can also be seen. Figures EJg) , (h), and (i) are those for 
N = 15212 at points 6 = 0.768, b = 0.769, and b = 0.77, 
respectively. We find from the figures that the transition 
point for the N — 15212 surface is close to b = 0.768 or 
6 = 0.769. 

The transition point moves left in the 6-axis with the 
increasing N as can be found in Figs. [Ja) and (b). The 
reason for this is because of the size effect. We also note 
that two distinct peaks can be observed in the distri- 
bution (or the histogram) h(X 2 ) of X 2 , which are not 
depicted here. On the smaller surfaces of N = 5672, 
N = 4842, and N = 3242, two distinct peaks can also 
be observed in h(X 2 ) at each transition point. 

Figure^a) is a normalized histogram h(S2) of Sz/Ng 
obtained at the transition point 6 = 0.772 on the surface 
of size N = 7292. We observed two clear peaks in the 
/i(S l 2) shown in Fig. 01a). Note that two distinct states, 
which differ in size, were observed at the point 6 = 0.772 
as shown in Fig. [2Jb). The variation of S2 against MCS 
is plotted in Fig. Efb). FiguresHc) and^d) show h(S 2 ) 
obtained at 6 = 0.77 on the N= 10242 surface, and Figs. 
H» andlHf) show those of N = 15212 at 6 = 0.768. We 
found two distinct peaks in h(S2) shown in Figs. EUc) and 



EJe), just like for N = 7292 in Fig. Eta). We found two 
peaks also in h(S2) of TV = 15212 at 6 = 0.769, which was 
not presented in a figure. The two peaks in h(S2) shown 
in Figs. EUa), CI{c) and|3fe) indicate that the model un- 
dergoes a first-order phase transition. Moreover, we find 
from these figures that the discontinuous nature of the 
transition is visible only on the surfaces of size N > 7292. 
We obtained h(S2) on N < 5762 surfaces and confirmed 
that there are no clear two-peaks in h(S2) in contrast to 
those shown in Figs. GIa),|3{c), and|3Je). 

In order to obtain the Hausdorff dimensions H snw in 
the smooth phase close to the transition point and H cm 
in the crumpled phase close to the transition point, the 
mean value of X 2 is obtained by averaging X 2 over a 
small region at each peak of h(X 2 ): 28 < X 2 < 80 and 
118<X 2 < 165 at 6 = 0.769 on the surface of TV = 15212, 
20<X 2 <55 and 78<X 2 <110 at 6 = 0.77 on iV=10242, 
15 < X 2 < 45 and 53 < X 2 < 82 at 6 = 0.772 on iV = 7292, 
12 < X 2 < 38 and 40 < X 2 < 68 at 6 = 0.773 on N = 5762, 
and 12 < X 2 < 30 and 35 < X 2 < 55 at 6 = 0.772 on 
TV = 4842. 

Figure^Ja) shows log-log plots of X 2 against N, which 
are obtained by averaging X 2 in the regions listed as 
above. Error bars on the data are the standard devi- 
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FIG. 4: (a) Log-log plots of X 2 against TV obtained in the 
smooth phase and in the crumpled phase at the transition 
point of surfaces TV > 4842. (b) Linear-log plots of X 2 against 
TV. The error bars on the data are the standard deviations. 



ations. The straight lines are drawn in the figures by 
fitting the data X 2 to 



X 2 ~ N 2/H , 

and as a consequence we have 

H smo = 2.02 ±0.14, 
H au = 2.59 ± 0.57. 



(4) 



(5) 



It is remarkable that H cru is less than the physical bound 
if = 3, although the error is relatively large. We note also 
that -f/cru in Eq. © is comparable to the theoretical pre- 
diction H = 2.39(23) within the error, which corresponds 
to the scaling exponent v — 0.84±0.04 [2^ where v — 2 /H. 
On the contrary, H smo is almost identical to the topolog- 
ical dimension 2. This indicates that the surface can be 
viewed as a smooth surface in the smooth phase. 

The large error in H crn indicates that H C1U is not well 
defined. Hence, it is possible that X 2 is logarithmically 
divergent |2^, 1^3, EH|- To check the logarithmic diver- 
gence of X 2 in the crumpled phase and in the smooth 
phase, we plot X 2 against log TV in Fig. Efb). We im- 
mediately find from the figure that X 2 in the smooth 
phase does not scale according to X 2 = cq + c\ log TV, 
which is expected just in the limit b—>0. It is also found 
that X 2 does not scale according to the logarithmic di- 
vergence in the crumpled phase. In fact, a dimension- 
less quantity the residual sum of squares RSS, defined 
by RSS = J^[(data - fitting formula) /error] 2 , for X 2 in 
the crumpled phase is RSS = 0.403 which is obtained 
by the linear-log fit in Fig. Efb), and it is larger than 
RSS = 0.165 which is obtained by the log- log fit in Fig. 
EKa). Thus, the log- log fit is better than the linear- log 
fit for X 2 in the crumpled phase at the transition point. 
This allows us to conclude that H cm in Eq. is a 
well-defined value. We must note, however, that only 
large scale simulations can clarify whether H cru has a 
well-defined value and is less than the physical bound. 

In order to see the critical slowing down typical to 
phase transitions, we plot in Figs. 03a) and |Sfb) the 



autocorrelation coefficient A(X 2 ) of X 2 defined by 
E^ 2 (r t )X 2 (r l+1 ) 



A{X Z ) 



n x 500, n=l,2, 



(6) 



Tj+l = Ti 

The horizontal axes in the figure represent 500 X n (n- 
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FIG. 5: A(X 2 ) obtained at the smooth phase, the crumpled 
phase, and the transition point on surfaces of (a) TV = 10242 
and (b)TV= 15212. 

1, 2, • • • )-MCS, which are a sampling-sweep between the 
samples X 2 {n) and X 2 (t;+i). The critical slowing down 
is clearly seen in the figure. The reason for this is because 
the volume of phase space (C R 3 ) of X at the transition 
point becomes larger than those at the crumpled phase 
and at the smooth phase. We find also an expected be- 
havior for A(X 2 ) such that A(X 2 ) depends on TV and 
A(X 2 ) becomes larger with increasing TV, which is not 
plotted in the figure. 

The specific heat Cs 2 defined by 



c S2 = b ^((s 2 -(s 2 )) 2 ) 



(7) 



can reflect the phase transition. Figure EH a) shows Cs 2 
which were obtained on surfaces of size TV = 15212, TV = 
10242, and TV = 2562 respectively. Sharp peaks of Cs 2 
in Fig. HJa) indicate a discontinuous transition. The 
peaks are located at b = 0.768, b = 0.77, b — 0.776 on 
surfaces of size TV = 15212, TV = 10242, and TV = 2562 
respectively. The curve of Cs 2 of TV = 2562 is relatively 
smooth, however, it is not so smooth for TV = 15212 or 
TV = 10242. These irregular behaviors for Cs 2 can be 
seen more or less equally for TV = 7292, TV = 5762 or 
TV = 4842. On these surfaces, we can see two peaks in the 
histograms h(X 2 ). Thus, we understand that it is very 
hard to obtain Cg 2 smoothly from such surfaces whose 
size changes discontinuously. 

Figure Efb) is a log-log plot of the peak value Cg^ ax 
against TV including the results obtained on surfaces of 
TV = 3242, TV = 2562, TV = 1442, and TV = 812. The straight 
lines were drawn by fitting the largest four Cc? ax and the 



smallest four C!? ax to 



a 



s 2 



(8) 
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FIG. 6: (a) The specific heat Cs 2 against the bending rigidity 
b, and (b) the peak values C<^ ax against the number of vertices 
iV in a log-log scale. The error bars on the data are the 
statistical errors. The straight lines in (b) are drawn by fitting 
the largest four and the smallest four Cgf* to Eq. ©. The 
units of Cs 2 and b are (kT) 2 and kT respectively. 

where a is a critical exponent of the transition. Thus, we 
have 

(Ti = 0.93 ±0.13, (N > 5762), 

cr 2 = 0.45 ±0.14, (N < 4842). (9) 

(Ti =0.93(13) which indicates that the phase transition is 
of the first order. On the contrary, 02 =0.45(14) in Eq. 
(|9"|) implies that the model appears to undergo a contin- 
uous transition at N < 4842. Two different behaviors of 
C!? ax against ./V shown in Fig. E^b) are consistent with 
the fact that two distinct peaks in h(S2) can be observed 
only on larger (TV > 7292) surfaces. 

Figure E^a) is a snapshot of the N — 15212 surface 
in the crumpled phase at b = 0.769, and Fig. [7£b) is 
the one in the smooth phase at the same b. The mean 
square size is about X 2 = 54 and X 2 = 138 in Figs. EJa) 
and [7{b) , respectively. The sections for the surfaces in 
Figs. EJa) and [7|[b) are depicted in Figs. |7Jc) and[7|[d) 
respectively. Surfaces are rough in short scales even in 
the smooth phase shown in Figs. Elb), whereas they are 
smooth in the long range scale. Surfaces rough in short 
scales can also be seen deep in the smooth phase. There 
appears to be only a spherical mono-layer surface in the 
smooth phase, and there are no apparent oblong, linear 
or branched polymer surfaces for N < 15212 at least. 
On the contrary, surfaces in the crumpled phase at the 
transition point are not completely collapsed, and they 
appear not only crumpled but also smooth. Thus, we 
understand a reason why the Hausdorff dimension H CTU 
in Eq. JSJ is less than the physical bound. 

Finally, we plot Si/N against b in Figs. OUa) and|SIb) 
in order to check that the equilibrium configurations were 
obtained in the MC simulations. We see that the ex- 
pected relation Si/N = 3(N- 1)/(2N) ~ 3/2 holds in all 
cases in the figure. The deviations are very small. It was 
also confirmed that the relation Si/N ~ 3/2 is satisfied 
in all other cases, which were not presented in the figure. 




FIG. 7: Snapshots of surfaces at (a) the crumpled phase and 
at (b) the smooth phase, and (c) the section of the surface in 
(a), and (d) the section of the surface in (b). The snapshots 
were obtained at 6 = 0.769 on the surface of iV = 15212. The 
mean square size is about (a) X 2 = 54 and (b) X 2 = 138. 



IV. SUMMARY AND CONCLUSIONS 

To conclude, we have shown that the tethered sur- 
face model of Helfrich and Polyakov-Kleinert undergoes 
a first-order phase transition between the smooth and 
the crumpled phases by performing canonical MC simu- 
lations on spherical surfaces of size up to N— 15212. The 
surfaces were constructed by dividing the icosahedron. 
The discrete form of the Hamiltonian, a linear combina- 
tion of the Gaussian term and the bending energy term, 
is the one widely used in the numerical studies carried 
out so far. The first-order transition was observed on 
surfaces of size N > 7292. We have checked that the tran- 
sition appears as a second-order one on relatively small 
surfaces, and hence that our results are consistent with 
results reported in previous numerical studies. We also 
confirmed that the Hausdorff dimension H smo is close to 
the topological dimension H — 2 in the smooth phase at 
the transition point, and that 7J cru remains finite and is 
less than the physical bound H — 3 in the crumpled phase 
at the transition point. Consequently, it is possible to 
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FIG. 8: Si/N against b obtained on surfaces of (a) N = 15212, 
7V = 10242, Af = 7292, and (b) Ar = 5762, 7V = 4842, AT = 3252. 



consider that the model represents a smooth and crum- 
pled state in real physical membranes at the crumpled 
phase close to the transition point. Further numerical 



studies would clarify this problem. Since the first-order 
nature of the transition has been confirmed in the HPK 
model, it would be interesting to study the model in more 
detail with large scale simulations: the phase diagram in 
the tensionless model, the dependence of the transition 
on the topology of surfaces, and the phase diagram of the 
fluid surfaces. 

The simulations were done by using a Pentium-4 
(2.8GHz-3.2GHz) PC with an Intel Fortran Compiler for 
Linux and one for Windows. The total number of CPU 
time was about 1350 days. Snapshots of surfaces were 
generated with POV-Ray for Windows v3.5. 
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